Source of the materials: Biopython cookbook (adapted) Status: Draft
Swiss-Prot and ExPASy¶
Parsing Swiss-Prot files¶
Swiss-Prot (http://www.expasy.org/sprot) is a hand-curated database of protein sequences. Biopython can parse the “plain text” Swiss-Prot file format, which is still used for the UniProt Knowledgebase which combined Swiss-Prot, TrEMBL and PIR-PSD. We do not (yet) support the UniProtKB XML file format.
Parsing Swiss-Prot records¶
In Section [sec:SeqIO_ExPASy_and_SwissProt], we described how to
extract the sequence of a Swiss-Prot record as a SeqRecord object.
Alternatively, you can store the Swiss-Prot record in a
Bio.SwissProt.Record object, which in fact stores the complete
information contained in the Swiss-Prot record. In this section, we
describe how to extract Bio.SwissProt.Record objects from a
Swiss-Prot file.
To parse a Swiss-Prot record, we first get a handle to a Swiss-Prot record. There are several ways to do so, depending on where and how the Swiss-Prot record is stored:
- Open a Swiss-Prot file locally:
>>> handle = open("myswissprotfile.dat") - Open a gzipped Swiss-Prot file:
In [2]:
import gzip
handle = gzip.open("myswissprotfile.dat.gz")
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
<ipython-input-2-5b562ec6e633> in <module>()
1 import gzip
----> 2 handle = gzip.open("myswissprotfile.dat.gz")
/home/tiago_antao/miniconda/lib/python3.5/gzip.py in open(filename, mode, compresslevel, encoding, errors, newline)
51 gz_mode = mode.replace("t", "")
52 if isinstance(filename, (str, bytes)):
---> 53 binary_file = GzipFile(filename, gz_mode, compresslevel)
54 elif hasattr(filename, "read") or hasattr(filename, "write"):
55 binary_file = GzipFile(None, gz_mode, compresslevel, filename)
/home/tiago_antao/miniconda/lib/python3.5/gzip.py in __init__(self, filename, mode, compresslevel, fileobj, mtime)
161 mode += 'b'
162 if fileobj is None:
--> 163 fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
164 if filename is None:
165 filename = getattr(fileobj, 'name', '')
FileNotFoundError: [Errno 2] No such file or directory: 'myswissprotfile.dat.gz'
- Open a Swiss-Prot file over the internet:
In [5]:
import urllib.request
handle = urllib.request.urlopen("http://www.somelocation.org/data/someswissprotfile.dat")
---------------------------------------------------------------------------
gaierror Traceback (most recent call last)
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
1239 try:
-> 1240 h.request(req.get_method(), req.selector, req.data, headers)
1241 except OSError as err: # timeout error
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in request(self, method, url, body, headers)
1082 """Send a complete request to the server."""
-> 1083 self._send_request(method, url, body, headers)
1084
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in _send_request(self, method, url, body, headers)
1127 body = body.encode('iso-8859-1')
-> 1128 self.endheaders(body)
1129
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in endheaders(self, message_body)
1078 raise CannotSendHeader()
-> 1079 self._send_output(message_body)
1080
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in _send_output(self, message_body)
910
--> 911 self.send(msg)
912 if message_body is not None:
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in send(self, data)
853 if self.auto_open:
--> 854 self.connect()
855 else:
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in connect(self)
825 self.sock = self._create_connection(
--> 826 (self.host,self.port), self.timeout, self.source_address)
827 self.sock.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)
/home/tiago_antao/miniconda/lib/python3.5/socket.py in create_connection(address, timeout, source_address)
692 err = None
--> 693 for res in getaddrinfo(host, port, 0, SOCK_STREAM):
694 af, socktype, proto, canonname, sa = res
/home/tiago_antao/miniconda/lib/python3.5/socket.py in getaddrinfo(host, port, family, type, proto, flags)
731 addrlist = []
--> 732 for res in _socket.getaddrinfo(host, port, family, type, proto, flags):
733 af, socktype, proto, canonname, sa = res
gaierror: [Errno -2] Name or service not known
During handling of the above exception, another exception occurred:
URLError Traceback (most recent call last)
<ipython-input-5-19b02eddd780> in <module>()
1 import urllib.request
----> 2 handle = urllib.request.urlopen("http://www.somelocation.org/data/someswissprotfile.dat")
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
160 else:
161 opener = _opener
--> 162 return opener.open(url, data, timeout)
163
164 def install_opener(opener):
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
463 req = meth(req)
464
--> 465 response = self._open(req, data)
466
467 # post-process response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _open(self, req, data)
481 protocol = req.type
482 result = self._call_chain(self.handle_open, protocol, protocol +
--> 483 '_open', req)
484 if result:
485 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_open(self, req)
1266
1267 def http_open(self, req):
-> 1268 return self.do_open(http.client.HTTPConnection, req)
1269
1270 http_request = AbstractHTTPHandler.do_request_
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
1240 h.request(req.get_method(), req.selector, req.data, headers)
1241 except OSError as err: # timeout error
-> 1242 raise URLError(err)
1243 r = h.getresponse()
1244 except:
URLError: <urlopen error [Errno -2] Name or service not known>
- Open a Swiss-Prot file over the internet from the ExPASy database (see section [subsec:expasy_swissprot]):
In [6]:
from Bio import ExPASy
handle = ExPASy.get_sprot_raw(myaccessionnumber)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-6-c988fb329f11> in <module>()
1 from Bio import ExPASy
----> 2 handle = ExPASy.get_sprot_raw(myaccessionnumber)
NameError: name 'myaccessionnumber' is not defined
The key point is that for the parser, it doesn’t matter how the handle was created, as long as it points to data in the Swiss-Prot format.
We can use Bio.SeqIO as described in
Section [sec:SeqIO_ExPASy_and_SwissProt] to get file format agnostic
SeqRecord objects. Alternatively, we can use Bio.SwissProt get
Bio.SwissProt.Record objects, which are a much closer match to the
underlying file format.
To read one Swiss-Prot record from the handle, we use the function
read():
In [7]:
from Bio import SwissProt
record = SwissProt.read(handle)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-7-4150bec69677> in <module>()
1 from Bio import SwissProt
----> 2 record = SwissProt.read(handle)
NameError: name 'handle' is not defined
This function should be used if the handle points to exactly one
Swiss-Prot record. It raises a ValueError if no Swiss-Prot record
was found, and also if more than one record was found.
We can now print out some information about this record:
In [8]:
print(record.description)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-8-040616e07cd0> in <module>()
----> 1 print(record.description)
NameError: name 'record' is not defined
In [9]:
for ref in record.references:
print("authors:", ref.authors)
print("title:", ref.title)
File "<ipython-input-9-a58c3c1279e1>", line 2
print("authors:", ref.authors)
^
IndentationError: expected an indented block
In [10]:
print(record.organism_classification)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-10-e3473e90e26a> in <module>()
----> 1 print(record.organism_classification)
NameError: name 'record' is not defined
To parse a file that contains more than one Swiss-Prot record, we use
the parse function instead. This function allows us to iterate over
the records in the file.
For example, let’s parse the full Swiss-Prot database and collect all
the descriptions. You can download this from the ExPAYs FTP
site
as a single gzipped-file uniprot_sprot.dat.gz (about 300MB). This is
a compressed file containing a single file, uniprot_sprot.dat (over
1.5GB).
As described at the start of this section, you can use the Python
library gzip to open and uncompress a .gz file, like this:
In [12]:
import gzip
handle = gzip.open("data/uniprot_sprot.dat.gz")
However, uncompressing a large file takes time, and each time you open
the file for reading in this way, it has to be decompressed on the fly.
So, if you can spare the disk space you’ll save time in the long run if
you first decompress the file to disk, to get the uniprot_sprot.dat
file inside. Then you can open the file for reading as usual:
In [13]:
handle = open("data/uniprot_sprot.dat")
As of June 2009, the full Swiss-Prot database downloaded from ExPASy contained 468851 Swiss-Prot records. One concise way to build up a list of the record descriptions is with a list comprehension:
In [15]:
from Bio import SwissProt
handle = open("data/uniprot_sprot.dat")
descriptions = [record.description for record in SwissProt.parse(handle)]
len(descriptions)
Out[15]:
549832
In [16]:
descriptions[:5]
Out[16]:
['RecName: Full=Putative transcription factor 001R;',
'RecName: Full=Uncharacterized protein 002L;',
'RecName: Full=Uncharacterized protein 002R;',
'RecName: Full=Uncharacterized protein 003L;',
'RecName: Full=Uncharacterized protein 3R; Flags: Precursor;']
Or, using a for loop over the record iterator:
In [17]:
from Bio import SwissProt
descriptions = []
handle = open("data/uniprot_sprot.dat")
for record in SwissProt.parse(handle):
descriptions.append(record.description)
In [18]:
len(descriptions)
Out[18]:
549832
Because this is such a large input file, either way takes about eleven
minutes on my new desktop computer (using the uncompressed
uniprot_sprot.dat file as input).
It is equally easy to extract any kind of information you’d like from Swiss-Prot records. To see the members of a Swiss-Prot record, use
In [19]:
dir(record)
Out[19]:
['__class__',
'__delattr__',
'__dict__',
'__dir__',
'__doc__',
'__eq__',
'__format__',
'__ge__',
'__getattribute__',
'__gt__',
'__hash__',
'__init__',
'__le__',
'__lt__',
'__module__',
'__ne__',
'__new__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'accessions',
'annotation_update',
'comments',
'created',
'cross_references',
'data_class',
'description',
'entry_name',
'features',
'gene_name',
'host_organism',
'host_taxonomy_id',
'keywords',
'molecule_type',
'organelle',
'organism',
'organism_classification',
'references',
'seqinfo',
'sequence',
'sequence_length',
'sequence_update',
'taxonomy_id']
Parsing the Swiss-Prot keyword and category list¶
Swiss-Prot also distributes a file keywlist.txt, which lists the
keywords and categories used in Swiss-Prot. The file contains entries in
the following form:
ID 2Fe-2S.
AC KW-0001
DE Protein which contains at least one 2Fe-2S iron-sulfur cluster: 2 iron
DE atoms complexed to 2 inorganic sulfides and 4 sulfur atoms of
DE cysteines from the protein.
SY Fe2S2; [2Fe-2S] cluster; [Fe2S2] cluster; Fe2/S2 (inorganic) cluster;
SY Di-mu-sulfido-diiron; 2 iron, 2 sulfur cluster binding.
GO GO:0051537; 2 iron, 2 sulfur cluster binding
HI Ligand: Iron; Iron-sulfur; 2Fe-2S.
HI Ligand: Metal-binding; 2Fe-2S.
CA Ligand.
//
ID 3D-structure.
AC KW-0002
DE Protein, or part of a protein, whose three-dimensional structure has
DE been resolved experimentally (for example by X-ray crystallography or
DE NMR spectroscopy) and whose coordinates are available in the PDB
DE database. Can also be used for theoretical models.
HI Technical term: 3D-structure.
CA Technical term.
//
ID 3Fe-4S.
...
The entries in this file can be parsed by the parse function in the
Bio.SwissProt.KeyWList module. Each entry is then stored as a
Bio.SwissProt.KeyWList.Record, which is a Python dictionary.
In [20]:
from Bio.SwissProt import KeyWList
handle = open("data/keywlist.txt")
records = KeyWList.parse(handle)
for record in records:
print(record['ID'])
print(record['DE'])
2Fe-2S.
Protein which contains at least one 2Fe-2S iron-sulfur cluster: 2 iron atoms complexed to 2 inorganic sulfides and 4 sulfur atoms of cysteines from the protein.
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
<ipython-input-20-8c149ecad36c> in <module>()
3 records = KeyWList.parse(handle)
4 for record in records:
----> 5 print(record['ID'])
6 print(record['DE'])
KeyError: 'ID'
This prints
2Fe-2S.
Protein which contains at least one 2Fe-2S iron-sulfur cluster: 2 iron atoms
complexed to 2 inorganic sulfides and 4 sulfur atoms of cysteines from the
protein.
...
Parsing Prosite records¶
Prosite is a database containing protein domains, protein families,
functional sites, as well as the patterns and profiles to recognize
them. Prosite was developed in parallel with Swiss-Prot. In Biopython, a
Prosite record is represented by the Bio.ExPASy.Prosite.Record
class, whose members correspond to the different fields in a Prosite
record.
In general, a Prosite file can contain more than one Prosite records.
For example, the full set of Prosite records, which can be downloaded as
a single file (prosite.dat) from the ExPASy FTP
site, contains
2073 records (version 20.24 released on 4 December 2007). To parse such
a file, we again make use of an iterator:
In [ ]:
from Bio.ExPASy import Prosite
handle = open("myprositefile.dat")
records = Prosite.parse(handle)
We can now take the records one at a time and print out some information. For example, using the file containing the complete Prosite database, we’d find
In [ ]:
from Bio.ExPASy import Prosite
handle = open("prosite.dat")
records = Prosite.parse(handle)
record = next(records)
record.accession
In [ ]:
record.name
In [ ]:
record.pdoc
In [ ]:
record = next(records)
record.accession
In [ ]:
record.name
In [ ]:
record.pdoc
In [ ]:
record = next(records)
record.accession
In [ ]:
record.name
In [ ]:
record.pdoc
and so on. If you’re interested in how many Prosite records there are, you could use
In [ ]:
from Bio.ExPASy import Prosite
handle = open("prosite.dat")
records = Prosite.parse(handle)
n = 0
for record in records:
n += 1
In [ ]:
n
To read exactly one Prosite from the handle, you can use the read
function:
In [ ]:
from Bio.ExPASy import Prosite
handle = open("mysingleprositerecord.dat")
record = Prosite.read(handle)
This function raises a ValueError if no Prosite record is found, and also if more than one Prosite record is found.
Parsing Prosite documentation records¶
In the Prosite example above, the record.pdoc accession numbers
'PDOC00001', 'PDOC00004', 'PDOC00005' and so on refer to
Prosite documentation. The Prosite documentation records are available
from ExPASy as individual files, and as one file (prosite.doc)
containing all Prosite documentation records.
We use the parser in Bio.ExPASy.Prodoc to parse Prosite
documentation records. For example, to create a list of all accession
numbers of Prosite documentation record, you can use
In [ ]:
from Bio.ExPASy import Prodoc
handle = open("prosite.doc")
records = Prodoc.parse(handle)
accessions = [record.accession for record in records]
Again a read() function is provided to read exactly one Prosite
documentation record from the handle.
Parsing Enzyme records¶
ExPASy’s Enzyme database is a repository of information on enzyme nomenclature. A typical Enzyme record looks as follows:
ID 3.1.1.34
DE Lipoprotein lipase.
AN Clearing factor lipase.
AN Diacylglycerol lipase.
AN Diglyceride lipase.
CA Triacylglycerol + H(2)O = diacylglycerol + a carboxylate.
CC -!- Hydrolyzes triacylglycerols in chylomicrons and very low-density
CC lipoproteins (VLDL).
CC -!- Also hydrolyzes diacylglycerol.
PR PROSITE; PDOC00110;
DR P11151, LIPL_BOVIN ; P11153, LIPL_CAVPO ; P11602, LIPL_CHICK ;
DR P55031, LIPL_FELCA ; P06858, LIPL_HUMAN ; P11152, LIPL_MOUSE ;
DR O46647, LIPL_MUSVI ; P49060, LIPL_PAPAN ; P49923, LIPL_PIG ;
DR Q06000, LIPL_RAT ; Q29524, LIPL_SHEEP ;
//
In this example, the first line shows the EC (Enzyme Commission) number of lipoprotein lipase (second line). Alternative names of lipoprotein lipase are “clearing factor lipase”, “diacylglycerol lipase”, and “diglyceride lipase” (lines 3 through 5). The line starting with “CA” shows the catalytic activity of this enzyme. Comment lines start with “CC”. The “PR” line shows references to the Prosite Documentation records, and the “DR” lines show references to Swiss-Prot records. Not of these entries are necessarily present in an Enzyme record.
In Biopython, an Enzyme record is represented by the
Bio.ExPASy.Enzyme.Record class. This record derives from a Python
dictionary and has keys corresponding to the two-letter codes used in
Enzyme files. To read an Enzyme file containing one Enzyme record, use
the read function in Bio.ExPASy.Enzyme:
In [ ]:
from Bio.ExPASy import Enzyme
with open("data/lipoprotein.txt") as handle:
record = Enzyme.read(handle)
In [ ]:
record["ID"]
In [ ]:
record["DE"]
In [ ]:
record["AN"]
In [ ]:
record["CA"]
In [ ]:
record["PR"]
In [ ]:
record["CC"]
In [ ]:
record["DR"]
The read function raises a ValueError if no Enzyme record is found,
and also if more than one Enzyme record is found.
The full set of Enzyme records can be downloaded as a single file
(enzyme.dat) from the ExPASy FTP
site, containing
4877 records (release of 3 March 2009). To parse such a file containing
multiple Enzyme records, use the parse function in
Bio.ExPASy.Enzyme to obtain an iterator:
In [ ]:
from Bio.ExPASy import Enzyme
handle = open("enzyme.dat")
records = Enzyme.parse(handle)
We can now iterate over the records one at a time. For example, we can make a list of all EC numbers for which an Enzyme record is available:
In [ ]:
ecnumbers = [record["ID"] for record in records]
Accessing the ExPASy server¶
Swiss-Prot, Prosite, and Prosite documentation records can be downloaded from the ExPASy web server at http://www.expasy.org. Six kinds of queries are available from ExPASy:
- get_prodoc_entry
- To download a Prosite documentation record in HTML format
- get_prosite_entry
- To download a Prosite record in HTML format
- get_prosite_raw
- To download a Prosite or Prosite documentation record in raw format
- get_sprot_raw
- To download a Swiss-Prot record in raw format
- sprot_search_ful
- To search for a Swiss-Prot record
- sprot_search_de
- To search for a Swiss-Prot record
To access this web server from a Python script, we use the
Bio.ExPASy module.
Retrieving a Swiss-Prot record¶
Let’s say we are looking at chalcone synthases for Orchids (see section [sec:orchids] for some justification for looking for interesting things about orchids). Chalcone synthase is involved in flavanoid biosynthesis in plants, and flavanoids make lots of cool things like pigment colors and UV protectants.
If you do a search on Swiss-Prot, you can find three orchid proteins for Chalcone Synthase, id numbers O23729, O23730, O23731. Now, let’s write a script which grabs these, and parses out some interesting information.
First, we grab the records, using the get_sprot_raw() function of
Bio.ExPASy. This function is very nice since you can feed it an id
and get back a handle to a raw text record (no HTML to mess with!). We
can the use Bio.SwissProt.read to pull out the Swiss-Prot record, or
Bio.SeqIO.read to get a SeqRecord. The following code accomplishes
what I just wrote:
In [21]:
from Bio import ExPASy
from Bio import SwissProt
In [22]:
accessions = ["O23729", "O23730", "O23731"]
records = []
In [23]:
for accession in accessions:
handle = ExPASy.get_sprot_raw(accession)
record = SwissProt.read(handle)
records.append(record)
---------------------------------------------------------------------------
IncompleteRead Traceback (most recent call last)
<ipython-input-23-37f8f432c5c4> in <module>()
1 for accession in accessions:
2 handle = ExPASy.get_sprot_raw(accession)
----> 3 record = SwissProt.read(handle)
4 records.append(record)
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/SwissProt/__init__.py in read(handle)
132 raise ValueError("No SwissProt record found")
133 # We should have reached the end of the record by now
--> 134 remainder = handle.read()
135 if remainder:
136 raise ValueError("More than one SwissProt record found")
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in read(self, amt)
444 else:
445 try:
--> 446 s = self._safe_read(self.length)
447 except IncompleteRead:
448 self._close_conn()
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in _safe_read(self, amt)
592 chunk = self.fp.read(min(amt, MAXAMOUNT))
593 if not chunk:
--> 594 raise IncompleteRead(b''.join(s), amt)
595 s.append(chunk)
596 amt -= len(chunk)
IncompleteRead: IncompleteRead(0 bytes read, 3386 more expected)
If the accession number you provided to ExPASy.get_sprot_raw does
not exist, then SwissProt.read(handle) will raise a ValueError.
You can catch ValueException exceptions to detect invalid accession
numbers:
In [24]:
for accession in accessions:
handle = ExPASy.get_sprot_raw(accession)
try:
record = SwissProt.read(handle)
except ValueException:
print("WARNING: Accession %s not found" % accession)
records.append(record)
---------------------------------------------------------------------------
IncompleteRead Traceback (most recent call last)
<ipython-input-24-e050fd24d4b9> in <module>()
3 try:
----> 4 record = SwissProt.read(handle)
5 except ValueException:
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/SwissProt/__init__.py in read(handle)
133 # We should have reached the end of the record by now
--> 134 remainder = handle.read()
135 if remainder:
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in read(self, amt)
445 try:
--> 446 s = self._safe_read(self.length)
447 except IncompleteRead:
/home/tiago_antao/miniconda/lib/python3.5/http/client.py in _safe_read(self, amt)
593 if not chunk:
--> 594 raise IncompleteRead(b''.join(s), amt)
595 s.append(chunk)
IncompleteRead: IncompleteRead(0 bytes read, 3386 more expected)
During handling of the above exception, another exception occurred:
NameError Traceback (most recent call last)
<ipython-input-24-e050fd24d4b9> in <module>()
3 try:
4 record = SwissProt.read(handle)
----> 5 except ValueException:
6 print("WARNING: Accession %s not found" % accession)
7 records.append(record)
NameError: name 'ValueException' is not defined
Searching Swiss-Prot¶
Now, you may remark that I knew the records’ accession numbers
beforehand. Indeed, get_sprot_raw() needs either the entry name or
an accession number. When you don’t have them handy, you can use one of
the sprot_search_de() or sprot_search_ful() functions.
sprot_search_de() searches in the ID, DE, GN, OS and OG lines;
sprot_search_ful() searches in (nearly) all the fields. They are
detailed on http://www.expasy.org/cgi-bin/sprot-search-de and
http://www.expasy.org/cgi-bin/sprot-search-ful respectively. Note that
they don’t search in TrEMBL by default (argument trembl). Note also
that they return HTML pages; however, accession numbers are quite easily
extractable:
In [25]:
from Bio import ExPASy
import re
In [26]:
handle = ExPASy.sprot_search_de("Orchid Chalcone Synthase")
# or:
# handle = ExPASy.sprot_search_ful("Orchid and {Chalcone Synthase}")
html_results = handle.read()
if "Number of sequences found" in html_results:
ids = re.findall(r'HREF="/uniprot/(\w+)"', html_results)
else:
ids = re.findall(r'href="/cgi-bin/niceprot\.pl\?(\w+)"', html_results)
---------------------------------------------------------------------------
HTTPError Traceback (most recent call last)
<ipython-input-26-ee45c5db7457> in <module>()
----> 1 handle = ExPASy.sprot_search_de("Orchid Chalcone Synthase")
2 # or:
3 # handle = ExPASy.sprot_search_ful("Orchid and {Chalcone Synthase}")
4 html_results = handle.read()
5 if "Number of sequences found" in html_results:
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/ExPASy/__init__.py in sprot_search_de(text, swissprot, trembl, cgi)
111 options = _urlencode(variables)
112 fullcgi = "%s?%s" % (cgi, options)
--> 113 handle = _urlopen(fullcgi)
114 return handle
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
160 else:
161 opener = _opener
--> 162 return opener.open(url, data, timeout)
163
164 def install_opener(opener):
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
507 if http_err:
508 args = (dict, 'default', 'http_error_default') + orig_args
--> 509 return self._call_chain(*args)
510
511 # XXX probably also want an abstract factory that knows when it makes
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_default(self, req, fp, code, msg, hdrs)
587 class HTTPDefaultErrorHandler(BaseHandler):
588 def http_error_default(self, req, fp, code, msg, hdrs):
--> 589 raise HTTPError(req.full_url, code, msg, hdrs, fp)
590
591 class HTTPRedirectHandler(BaseHandler):
HTTPError: HTTP Error 301: Moved Permanently
Retrieving Prosite and Prosite documentation records¶
Prosite and Prosite documentation records can be retrieved either in HTML format, or in raw format. To parse Prosite and Prosite documentation records with Biopython, you should retrieve the records in raw format. For other purposes, however, you may be interested in these records in HTML format.
To retrieve a Prosite or Prosite documentation record in raw format, use
get_prosite_raw(). For example, to download a Prosite record and
print it out in raw text format, use
In [27]:
from Bio import ExPASy
handle = ExPASy.get_prosite_raw('PS00001')
text = handle.read()
print(text)
---------------------------------------------------------------------------
HTTPError Traceback (most recent call last)
<ipython-input-27-8c96cea49d6d> in <module>()
1 from Bio import ExPASy
----> 2 handle = ExPASy.get_prosite_raw('PS00001')
3 text = handle.read()
4 print(text)
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/ExPASy/__init__.py in get_prosite_raw(id, cgi)
62 For a non-existing key, ExPASy returns nothing.
63 """
---> 64 return _urlopen("%s?%s" % (cgi, id))
65
66
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
160 else:
161 opener = _opener
--> 162 return opener.open(url, data, timeout)
163
164 def install_opener(opener):
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
507 if http_err:
508 args = (dict, 'default', 'http_error_default') + orig_args
--> 509 return self._call_chain(*args)
510
511 # XXX probably also want an abstract factory that knows when it makes
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_default(self, req, fp, code, msg, hdrs)
587 class HTTPDefaultErrorHandler(BaseHandler):
588 def http_error_default(self, req, fp, code, msg, hdrs):
--> 589 raise HTTPError(req.full_url, code, msg, hdrs, fp)
590
591 class HTTPRedirectHandler(BaseHandler):
HTTPError: HTTP Error 301: Moved Permanently
To retrieve a Prosite record and parse it into a Bio.Prosite.Record
object, use
In [28]:
from Bio import ExPASy
from Bio import Prosite
handle = ExPASy.get_prosite_raw('PS00001')
record = Prosite.read(handle)
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
<ipython-input-28-f4e345e9eb60> in <module>()
1 from Bio import ExPASy
----> 2 from Bio import Prosite
3 handle = ExPASy.get_prosite_raw('PS00001')
4 record = Prosite.read(handle)
ImportError: cannot import name 'Prosite'
The same function can be used to retrieve a Prosite documentation record
and parse it into a Bio.ExPASy.Prodoc.Record object:
In [29]:
from Bio import ExPASy
from Bio.ExPASy import Prodoc
handle = ExPASy.get_prosite_raw('PDOC00001')
record = Prodoc.read(handle)
---------------------------------------------------------------------------
HTTPError Traceback (most recent call last)
<ipython-input-29-7482db7f5b87> in <module>()
1 from Bio import ExPASy
2 from Bio.ExPASy import Prodoc
----> 3 handle = ExPASy.get_prosite_raw('PDOC00001')
4 record = Prodoc.read(handle)
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/ExPASy/__init__.py in get_prosite_raw(id, cgi)
62 For a non-existing key, ExPASy returns nothing.
63 """
---> 64 return _urlopen("%s?%s" % (cgi, id))
65
66
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
160 else:
161 opener = _opener
--> 162 return opener.open(url, data, timeout)
163
164 def install_opener(opener):
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
501 http_err = 0
502 args = (dict, proto, meth_name) + args
--> 503 result = self._call_chain(*args)
504 if result:
505 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_302(self, req, fp, code, msg, headers)
684 fp.close()
685
--> 686 return self.parent.open(new, timeout=req.timeout)
687
688 http_error_301 = http_error_303 = http_error_307 = http_error_302
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in open(self, fullurl, data, timeout)
469 for processor in self.process_response.get(protocol, []):
470 meth = getattr(processor, meth_name)
--> 471 response = meth(req, response)
472
473 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_response(self, request, response)
579 if not (200 <= code < 300):
580 response = self.parent.error(
--> 581 'http', request, response, code, msg, hdrs)
582
583 return response
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in error(self, proto, *args)
507 if http_err:
508 args = (dict, 'default', 'http_error_default') + orig_args
--> 509 return self._call_chain(*args)
510
511 # XXX probably also want an abstract factory that knows when it makes
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
441 for handler in handlers:
442 func = getattr(handler, meth_name)
--> 443 result = func(*args)
444 if result is not None:
445 return result
/home/tiago_antao/miniconda/lib/python3.5/urllib/request.py in http_error_default(self, req, fp, code, msg, hdrs)
587 class HTTPDefaultErrorHandler(BaseHandler):
588 def http_error_default(self, req, fp, code, msg, hdrs):
--> 589 raise HTTPError(req.full_url, code, msg, hdrs, fp)
590
591 class HTTPRedirectHandler(BaseHandler):
HTTPError: HTTP Error 301: Moved Permanently
For non-existing accession numbers, ExPASy.get_prosite_raw returns a
handle to an emptry string. When faced with an empty string,
Prosite.read and Prodoc.read will raise a ValueError. You can
catch these exceptions to detect invalid accession numbers.
The functions get_prosite_entry() and get_prodoc_entry() are
used to download Prosite and Prosite documentation records in HTML
format. To create a web page showing one Prosite record, you can use
In [ ]:
from Bio import ExPASy
handle = ExPASy.get_prosite_entry('PS00001')
html = handle.read()
output = open("myprositerecord.html", "w")
output.write(html)
output.close()
and similarly for a Prosite documentation record:
In [30]:
from Bio import ExPASy
handle = ExPASy.get_prodoc_entry('PDOC00001')
html = handle.read()
output = open("myprodocrecord.html", "w")
output.write(html)
output.close()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-30-1fead8290c3a> in <module>()
3 html = handle.read()
4 output = open("myprodocrecord.html", "w")
----> 5 output.write(html)
6 output.close()
TypeError: write() argument must be str, not bytes
For these functions, an invalid accession number returns an error message in HTML format.
Scanning the Prosite database¶
ScanProsite allows you to scan protein sequences online against the Prosite database by providing a UniProt or PDB sequence identifier or the sequence itself. For more information about ScanProsite, please see the ScanProsite documentation as well as the documentation for programmatic access of ScanProsite.
You can use Biopython’s Bio.ExPASy.ScanProsite module to scan the
Prosite database from Python. This module both helps you to access
ScanProsite programmatically, and to parse the results returned by
ScanProsite. To scan for Prosite patterns in the following protein
sequence:
MEHKEVVLLLLLFLKSGQGEPLDDYVNTQGASLFSVTKKQLGAGSIEECAAKCEEDEEFT
CRAFQYHSKEQQCVIMAENRKSSIIIRMRDVVLFEKKVYLSECKTGNGKNYRGTMSKTKN
you can use the following code:
In [31]:
sequence = "MEHKEVVLLLLLFLKSGQGEPLDDYVNTQGASLFSVTKKQLGAGSIEECAAKCEEDEEFT
File "<ipython-input-31-a3137c4e2fe3>", line 1
sequence = "MEHKEVVLLLLLFLKSGQGEPLDDYVNTQGASLFSVTKKQLGAGSIEECAAKCEEDEEFT
^
SyntaxError: EOL while scanning string literal
In [32]:
from Bio.ExPASy import ScanProsite
handle = ScanProsite.scan(seq=sequence)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-32-fbf16249a192> in <module>()
1 from Bio.ExPASy import ScanProsite
----> 2 handle = ScanProsite.scan(seq=sequence)
NameError: name 'sequence' is not defined
By executing handle.read(), you can obtain the search results in raw
XML format. Instead, let’s use Bio.ExPASy.ScanProsite.read to parse
the raw XML into a Python object:
In [33]:
result = ScanProsite.read(handle)
type(result)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-33-8794b65cb5b5> in <module>()
----> 1 result = ScanProsite.read(handle)
2 type(result)
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/ExPASy/ScanProsite.py in read(handle)
65 saxparser = Parser()
66 saxparser.setContentHandler(content_handler)
---> 67 saxparser.parse(handle)
68 record = content_handler.record
69 return record
/home/tiago_antao/miniconda/lib/python3.5/xml/sax/expatreader.py in parse(self, source)
108 self.reset()
109 self._cont_handler.setDocumentLocator(ExpatLocator(self))
--> 110 xmlreader.IncrementalParser.parse(self, source)
111
112 def prepareParser(self, source):
/home/tiago_antao/miniconda/lib/python3.5/xml/sax/xmlreader.py in parse(self, source)
125 self.feed(buffer)
126 buffer = file.read(self._bufsize)
--> 127 self.close()
128
129 def feed(self, data):
/home/tiago_antao/miniconda/lib/python3.5/xml/sax/expatreader.py in close(self)
220 return
221 try:
--> 222 self.feed("", isFinal = 1)
223 self._cont_handler.endDocument()
224 self._parsing = 0
/home/tiago_antao/miniconda/lib/python3.5/site-packages/Bio/ExPASy/ScanProsite.py in feed(self, data, isFinal)
85 # fed to the parser.
86 if self.firsttime:
---> 87 if data[:5].decode('utf-8') != "<?xml":
88 raise ValueError(data)
89 self.firsttime = False
AttributeError: 'str' object has no attribute 'decode'
A Bio.ExPASy.ScanProsite.Record object is derived from a list, with
each element in the list storing one ScanProsite hit. This object also
stores the number of hits, as well as the number of search sequences, as
returned by ScanProsite. This ScanProsite search resulted in six hits:
In [ ]:
result.n_seq
In [ ]:
result.n_match
In [ ]:
len(result)
In [ ]:
result[0]
In [ ]:
result[1]
In [ ]:
result[2]
In [ ]:
result[3]
In [ ]:
result[4]
In [ ]:
result[5]
Other ScanProsite parameters can be passed as keyword arguments; see the
documentation for programmatic access of
ScanProsite
for more information. As an example, passing lowscore=1 to include
matches with low level scores lets use find one additional hit:
In [ ]:
handle = ScanProsite.scan(seq=sequence, lowscore=1)
result = ScanProsite.read(handle)
result.n_match